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INTRODUCTION 



The research program carried out under this research grant (WP-OIO67) was 
a two-year continuation of the studies initiated under a previous U.S. Public 
Health Service grant (WP-0079^). The final report describing the previous work 
was prepared by Noble (I966). The previous report shows the failure of Ekman- 
type models in synthesizing current meter records from buoy stations when using 
the wind measured at the station as an input stress. The previous study further 
suggests that geostrophic forces defined by the temperature-dependent density 
structure of Lake Michigan may be the dominant factor in determining the three- 
dimensional circulation patterns of the lake. 

Four major field experiments were carried out (in cooperation with other 
projects and individuals) under this grant. The two most significant operations 
were thermal mapping flights carried out by the Ant i- Submarine Warfare Environ- 
mental Prediction Services (ASWEPS) of the U.S. Naval Oceanographic Office. The 
results of the l8-21 October I966 flight program of the ASWEPS Super Constella- 
tion^ EL COYOTE (showing the existence of a major circulation feature^ the 
Coyote Current) have been described in detail in a report by Nobel^ £t al. (1967). 
The results of the second flighty on 25 October 1967^ are compared with those 
of the 1966 flight in a paper in preparation by V. E. Noble and J. C. Wilkerson. 

A third major field experiment was carried out in August 1967j> in coopera- 
tion with Dr. G. E. Birchfield of Northwestern University. A summary of the 
results of this study was presented by Birchfield and Noble. -^ In this study^ 
two current meter stations were established in Lake Michigan for a two-week 
period^ and comparisons were made between the current meter records, drogue 
measurements, and dynamic height current calculations. Spectral analyses and 
inertial period averages were also computed from the current meter records. 

The fourth field experiment consisted of a series of measurements of the 
temperature, eddy viscosity, eddy diffusivity, and current structures in the 
nearshore edge of the lake during the spring warming period of I968. These 
measurements were carried out by J. C. Huang. The results will be used as part 
of the basic input data and parameters for his theoretical investigations of 
lake circulation, and will be compared with the results obtained in connection 
with another project by Noble and Anderson (1968) during the spring of I967. 

In addition to the experimental program, two theoretical studies have been 
carried out in connection with this grant by students of the Department of 



Statistical analysis of currents at two nearby stations in Lake Michigan, summer 
1967.. G. E. Birchfield and V. E. Noble. Paper presented at 11th Conf . on 
Great Lakes Res., Milwaukee, I968. 



Meteorology and Oceanography of The University of Michigan. The first was the 
output from a Special Problems' course in which James H. Saylor calculated 
numerical results, based on Lake Michigan data, from Stern's (1965) model for 
the interaction of a uniform wind stress with a geostrophic vortex. The second 
is a Ph.D. dissertation being written by Joseph C. Huang. Mr. Huang is doing 
a theoretical analysis of a simplified Lake Michigan model using the early 
spring heat input to the lake as the prime current- gene rating force. Early re- 
sults from this model (Huang's dissertation is not expected to be completed 
until the summer of I969) indicate that the geostrophic currents resulting from 
the spring temperature structure of the lake may be as great as previous compu- 
tations for wind-driven currents. Those results correspond with observations, 
and with the concepts within Stern's (I965) model. 

In the proposal for this grant, it was anticipated that the data from a 
special triangular array of current meter stations set in Lake Huron during the 
winter of I965-66 would provide a basis for investigation of the Eulerian- 
Lagrangian transformation inherent in the extrapolation of current meter data, 
and of correlations between data from pairs of current meters. Instrumentation 
problems caused a lack of data from this source, and therefore it was necessary 
to rely upon measurements made under this grant for subsequent analysis and in- 
terpretation. 
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I. CHAEACTERIZATION OF THE CIRCULATION DYNAMICS OF LAKE MICHIGAN 

Vincent E. Noble 



The results of previous investigations carried out in connection with this 
program tend to indicate that the circulation dynamics of Lake Michigan and, by 
inference, the other Great Lakes, are dominated by geostrophic forces rather 
than the wind stresses during the summer season. Analysis of temperature records 
from the Great Lakes-Illinois River Basins Project thermographs installed at 
buoy stations in Lake Michigan during the winter season of 1965-6^ indicated 
that the prevailing westerly wind drift caused an eastward surface current on 
Lake Michigan, a sinking on the east side of the lake basin, a westward bottom 
current, and an upwelling on the west side of the basin (Heap and Noble I966) . 
In the spring, the "thermal bar" (Rodgers I965; Noble and Anderson I968), which 
is a consequence of the spring heating of the lake, defines a general circula- 
tion around the edge of the lake basin. During the onset of the summer season, 
with the development of the thermocline, the wind stress-geostrophic stress 
interactions become extremely complex. 

Critical testing of theoretical models depends heavily upon our ability to 
make adequate density measurements in both time and space. Analysis of 
data from single current meter stations has failed to demonstrate a direct rela- 
tionship between wind and current. Not enough is yet known about the spatial 
structure of currents to interpret the relationship between data from current 
meter stations separated by distances of 5 to 10 miles. Surface temperature 
patterns obtained by synoptic observations from aircraft imply certain types of 
circulation features when interpreted with respect to dynamic height current 
computations made from BT transects. However, this interpretation is subject 
to errors inasmuch as the dynamic height computations depend on the essential 
assumption that the current field is defined by the density (temperature) field, 
and therefore it is only self-consistent that the dynamic height currents should 
show a strong correlation with the surface temperature structure. To the extent 
that dynamic height computations show a good qualitative correlation with cur- 
rents measured with drogues and current meters, the fundamental assumption of 
the dominance of the geostrophic forces is valid. The question yet to be re- 
solved is the interaction between the wind and geostrophic forces. The persis- 
tence of small details in the surface temperature structure of the lake for 
periods of three weeks (Heap and Noble I966; Noble I967) give credence to the 
postulate of dominance of the lake circulation dynamics by geostrophic forces. 

Figure 1-1 shows isotherms contoured from 5-mile averages of surface 
temperature along 25 east-west flight lines spaced at 10-mile intervals from 
the south end of Lake Michigan on 18 October I966. Subsequent flight operations 
showed that the general temperature pattern persisted for at least three subse- 
quent days. Figure 1-2 shows the dynamic height currents computed for the 0-10 m 
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Figure 1-1. Sairface temperature contours from Airborne Radiation 
Thermometer flight of 18 October I966. 
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Figure 1-2. Surface temperatures and 0-10 m dynamic height currents 
from BT transects taken by R/v INLAND SEAS, October I966. 



layer, compared with surface temperature readings taken from BT transects along 
the 45 '^ 57, 2* N flight line on 18, 19, and 20 October I966 by R/v INIAp SEAS. 
These dynamic height currents indicated a strong northward current parallel to 
the strong temperature gradient on the west side of the northern half of the 
lake. This circulation feature was also indicated by a distinct visual water 
mass separation along the gradient. The water was warm and green to the east 
of the gradient, and cold and blue to the west. The circulation feature ^as 
named the Coyote Current in honor of the NAVOCEANO research aircraft, the Super 
Constellation EL COYOTE. 

If the circulation of Lake Michigan is dominated by geostrophic (temperature- 
dependent density) forces, it would be expeated that the Coyote Current would 
appear each fall in the northern end of Lake Michigan. Therefore another flight 
was scheduled during the fall of I967 to see if the current did recur. 

The early fall period of I967 was much colder than that of 1966, and the 
temperature mapping flight could not be scheduled until 25 October. Therefore 
the 1967 flight was equivalent to being two weeks later in the season than that 
of 1966. The surface temperature contours from the I967 Coyote flight are shown 
in Figure I-5. Due to heavy weather, the INLAKD SEAS could not make a BT transeci 
on the day of the flight. Figure l-k shows 0-10 m dynamic height currents and 
the surface temperature profile from the INLAM) SEAS BT transect at ky^^J.^^'N 
on 26 October I967. 

The 1967 temperature gradient was stronger (of the order of ^''C). The 
temperature gradient turned offshore near Sheboygan, Wis., in I966, and more 
nearly paralleled the shoreline in I967. Figure 1-5 shows the temperature struc- 
ture of a vertical section of the lake as determined by the BT transects of 
1966 and 1967- Because the development of the surface temperature structure 
was not followed continuously through the fall and early winter periods of 1966 
and 1967^ it is not possible to make definitive judgments based on the compari- 
son of measurements made on 18 October I966 and 25 October I967. However, the 
appearance of the strong temperature gradient feature covering a large portion 
of the north- south extent of the lake basin, accompanied by a water mass separa- 
tion as evidenced by color changes visible from the aircraft, gives a strong 
indication that a similar phenomenon was observed in both years. 

Additional small-scale experiments were planned and executed to provide 
experimental foundations for testing and evaluation of theoretical models that 
would lead to an understanding of the circulation dynamics of the lake. These 
experiments were continued observations of the current regime associated with 
the early spring temperature structure, and a 2-week period of measurement using 
two current meter stations and comparison current drogue and BT observations. 




Figure 1-5 • S-urface temperature contoiirs from Airborne Radiation 
Thermometer flight of 25 October I967. 
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Figure 1-k. Surface temperatures and 0-10 m dynamic height currents 
from BT transects taken by r/v INLAND SEAS, 26 October 196?. 
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Figure 1-5. Temperature structure of a vertical cross-section of Lake 
Michigan corresponding to Airborne Radiation Thermometer flight tracks, 
(a) 18 October 1966. (b) 20 October 1966. (c) 26 October I967. 
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II. CURRENT METER STUDIES OFF SHEBOYGAN, WISCONSIN, AUGUST I967 

Vincent E, Noble 



During the period of I6 August through 2k August IS&J , two current meter 
buoy stations were established in 280 ft of water about 8 miles offshore from 
Sheboygan, Wis. The station locations and relative bottom topography in that 
area of the lake are shown in Figure 2-1. The buoy stations were set a half 
mile apart in roughly an east-west direction with current meters set nominally 
at the surface and at 10-m depths. The actual current speed rotors were at 2 
and 12 m, respectively. Both current meter stations were obtained on loan from 
the Federal Water Pollution Control Administration, Great Lakes- Illinois River 
Basins Project. The experiment was designed to compare current drogue obser- 
vations and dynamic height current computations with current meter readings, and 
to provide an insight into the distance through which current meter readings 
can be extrapolated and into the Eulerian-Lagrangian transformation represented 
by current drogue and current meter measurements. 

Station I was set at l^i-i^O EST on 16 August I967 at a position of kykh.VE 
and 87''32.0'W. The station depth was 2&\- ft. Station II was set at I530 on 
16 August 1967 at a position hyhk.2''^ and 87^52. 6'W. The water depth at sta- 
tion II was 280 ft. The station position was one-half mile at bearing 287*^ from 
station I. The movie film on which the data from the surface meter at station 
I were recorded was prematurely exposed, therefore, no useful data were obtained 
for the first ^1-5 hours after the station was installed. The other three current 
meters operated satisfactorily throughout the full period. 

The buoy station configuration was similar to that used by the Great Lakes- 
Illinois River Basins Project. In this instance, however, the subsurface buoy 
(which is used to provide tension on the meter string) was set so that it was 
just awash at the surface of the water. The thermocline was rather shallow in 
the region of measurement so that the 12-m current rotor was just about at the 
thermocline, generally being in the metalimnion. The weather during the period 
of measurement was generally unfavorable so that far less comparison data were 
obtained than had been desired. On I7 August only BT's were taken for dynamic 
height current computations. On the l8th only drogue observations were made. 
On the 19th the ship was weathered-in at the dock. On the 20th both drogue and 
BT observations were made. On the 21st, BT*s only. The ship was laid in on the 
22d, with only BT*s again on the 25d. The stations were recovered on the 2^th. 
The winds were generally southwest to west on the l6th and 17th, with a frontal 
passage on the l8th producing northwest winds through the 19th. The wind went 
back to southwest to west on the 20th with a second cold front passage with 
north to northeast winds on the 21st remaining northeast through the 22d and 
going southeast on the 25d. 
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Figure 2-1. Location of current meter stations, showing 
contours of bottom topography. 
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The full statistical analyses of the current meter data have been discussed 
by Birchfield and Noble, The statistical analyses^ performed by Dr. G. E* 
Birchfield, Department of Engineering Sciences, Northwestern University, include 
determination of the mean transport current and a measure of the dispersion of 
the instantaneous current vectors during full and quarter- inert ial periods (17-I/5 
hours, and ^-l/5 hours, respectively). In this report, the quarter-period means 
are compared with drogue measurements, dynamic height current computations, and 
wind records from Milwaukee and Muskegon. The wind records from the land stations 
are used to document the general wind regime and significant frontal passages. 
Figure 2-2 compares the wind records from the land stations with the quarter- 
period averages from each of the current meters. 

The quarter-period current meter averages discussed in this section of the 
report are the result of a first analysis carried out by Dr. Birchfield. The 
quarter-period averages were computed by taking the times that the current meters 
were set as t = 0. Therefore, the 4-l/5-hour averages for stations I and II 
are displaced in time by 50 minutes (the elapsed time between the setting of 
stations I and II) . The current meter averages are being calculated for coin- 
cident periods, and a complete report of the statistical analysis of the current 
meter data and drogue comparisons is being prepared by Birchfield and Noble. 
The complete report will discuss the current meter averages and vector disper- 
sions for the coincident periods, and will contain spectral analyses of the 
current meter data. 

The winds at Milwaukee (MKE) were S from I5 August through 0000 I8 August. 
A frontal passage occurred around O5OO on I8 August, with winds shifting to N 
at 15 kts (ships on the lake reported 18-25 kts, N) . The MKE winds diminished 
on the 19th and went light, SW on the 20th (a ship reported JO i^"ts at 190° at 
2100 on 20 August). MKE winds were W-SW, 5-10 kts through 1200 21 August, with 
a strong frontal passage around I5OO on 21 August (ships reported 15-55 kts at 
060°-020° behind the front during the period 0900-2100 on 21 August). MKE 
winds went east, decreasing through I80O on 2k August. (Ship reports gave 21 kts 
at 030° for 2100 on 22 August, and 09 kts at 010° for 2100 on 25 August.) The 
ship reports used for this study were carferry reports taken approximately 50 
miles NE of the buoy stations. 

As shown in Figure 2-2, the quasi- steady currents between the fronts do not 
show a simple relation to the wind direction. The response of the current meter 
records to a frontal passage is the development of a rotational characteristic 
lasting for a period of I8 hours. (The rotational characteristic develops after 
a time lag of the order of 12 hours.) It is suggested that the complex relation . 
between the current vectors and wind field may be explained in Ekman transport 
and acceleration of geostrophic eddies as described by Stern ^s (I965) model. 



2 
See footnote 1, p. 1. 
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Figure 2-2 Comparison of Milwaukee (MKE) winds and 4-1/2 hour current vector 
averages for stations I and II at 0-10 m depths (15 August - 24 August 1967) 



Figures 2-5 through 2-7 show the results of the comparative measurements 
between the drogues, dynamic height current computations, and the quarter iner- 
tial period averages of the cu!i:rent meter records. The quarter period averages 
give the average current vector determined from a 4-l/5-hour interval. The 
dynamic height current computations yield currents between pairs of BT*s taken 
within a time interval of the order of l/2 hour within the comparison current 
meter period. The net drogue movements indicated by the dotted lines on the 
figures are used to determine a net transport vector. 

Figure 2-5 shows the relation between dynamic height currents and quarter 
period averages of the current meter records for 17 August 1967* For station 
I the averaging interval was O80O to 1220 hours and for station II, 085O to 
1510 on 17 August (period 5) The surface current data for station I were lost 
due to premature exposure of the film but records are available for 10 m at 
station I, and and 10 m at station II. The average currents from stations I 
and II are 2O-5I cm/sec. Dynamic height currents are computed from three BT 
casts taken between 0850 and 0952 on I7 August. The dynamic height computations 
give only the component of current speeds normal to the path between BT casts. 
For this comparison period the components determined by dynamic heights were 
k.9 and 5*6 cm/sec, respectively, for the average current in the to 10-m 
level between the three BT casts. 

The dynamic height currents were calculated by the method of Ayers (1957), 
and assume a reference depth of 70 m. The computed speed components were within 
20^ of those shown by the average values for the current meter stations. The 
dynamic height component computed for a position 2 miles east of the stations 
showed good agreement with the current meter average, but the component computed 
for a position I-I/2 miles west of the stations showed a direction I80** from 
that given by the station data. The ship winds at 0815 on I7 August were I9 
kts, with direction 195 °« 

Figure 2-k shows a comparison between drogue observations and quarter period 
current meter averages for 18 August 1967. The four current drogues were set 
between 0955 and IO55 and recovered between 1^^-15 and 15^1-1 on I8 August. A frontal 
passage had occurred Just before setting the drogues, and only two positions of 
the drogues were determined — the set point and the recovery point. Drogue speeds 
were generally to the north with a net transport velocity ranging from 5-7 to 
21.5 cm/sec. For station I a quarter period average encompassed the time be- 
tween 1000 and 1^20 on I8 August with surface current of 19*5 cm/sec and 10-m 
current of 55*6 cm/sec (period 11). The surface current speed was 5O.8 cm/sec 
with 10-m speed of 52.8 cm/sec at station II. 

The current drogues consisted of ^-ft square sheet metal current crosses 
suspended from modified net stake floats. The center of the current cross was 
at about the same depth as the current speed rotor for the surface meters. The 
drogue showed a good correlation with both speed and direction with the surface 
current meter of station II but very poor correlation with the direction of 
station I. The winds were generally north and at lUoO the ship wind was north 
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Figure 2-k. Drogue observations and quarter- period current meter 
averages, 18 August I968. The net drogue displacements and average 
transport velocities are shovra for a 5-hour interval. Wind N, 19- 
25 kts. 
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Figure 2-6. Dynamic height current components and corresponding 
i^-l/5 hour current meter averages, 21 August I967. Current speed 
in cm/sec. Wind S, 7 kts. 



21 



LEGEND 



CURRENT METER - Surface 
CURRENT METER - 10 meter 
-T' DROGUE 
— X DROGUE DISPLACEMENT 



a* DYNAMIC HEIGHTS 

O BT 



51.4 
/ 




25.9 21.3 



4.2 



\= 



43°46' 



O 



43044. 



33.7 



1/2 



• miles 



86° 34' 



43° 42' 



86°30' 



Figure 2-7. Dynamic height cairrent components and U-I/5 hour current 
meter averages, 23 August I967. The dynamic height coraponents were 
computed for a time at the end of averaging period 57 (solid vectors) 
and the beginning of period 58 (broken vectors). Wind SE, 15 kts. 
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19-25 kts. It is to be noted that all four drogues, and the surface current 
averaged at station II for the quarter period, were all to the north, that is, 
into the wind. 

Figures 2-5^ and b show a comparison between drogues, dynamic height cur- 
rents, and current meter averages for 20 August 1967* The five drogues were set 
on an east-west line through the current meter stations between 0912 and 0935* 
The drogues were fixed between 1525 and l4^0 and recovered between I655 sind 17^0. 

Figure 2- 5a shows the first half of the measurements between the set and 
the first fix of the drogues, and Figure 2-5b shows the second half between the 
fix and recovery. Figure 2- 5a gives the quarter period average currents for 
the time period between 09^0 and 1^00 for station I (period 22) and between 
1050 and 1^50 for station II (period 22). The dynamic height current component 
for a position 1 mile west of the current meter station shows good agreement 
(speed: 0.7 cm/sec) for the northerly component with the net transport of the 
drogue set at the same location (8.6 cm/sec to the north of east). The speed 
component (to the north at 0.1 cm/sec) computed for a position I-I/2 miles to 
the east of the stations is in the proper sense, but does not agree as well 
with its corresponding net drogue transport (9.8 cm/sec to the northeast). Of 
particular interest is the net south transport of the drogue set 2 miles to the 
west of the stations. There is good agreement between the drogues set at the 
current meter sites and the ^-l/j-hour meter averages. 

Figure 2- 5b compares the quarter period current meter records from 1^00 to 
1820 for station I (period 25) and 1^-50 to I9IO for station II (period 25) with 
the net transport currents given by the motion of the drogue between their fix 
and recovery points. These two direct current measurements are compared with 
the to 10-m current calculated by the dynamic height method. In this case, 
there is less correlation between the surface current meter averages, and be- 
tween these averages and the drogue transports and the dynamic height current 
components. (Only two dynamic height components are calculated for this inter- 
val, shown as 2.5 and 0.8 cm/sec to the east.) 

Figure 2-6 gives the results of the experiments on 21 August I967 comparing 
the 10-m dynamic height currents with the quarter period average from station I 
in the interval of 0720 to 11^0 (period 27) and the interval of 08IO to I25O 
at station II (period 27). Figure 2-7 gives the results for 25 August in which 
dynamic height currents for the to 10-m layer are compared with current meter 
averages for two adjacent U-hour periods. The BT casts for the computation of 
the dynamic height currents were taken at the end of period 57 in the averaging 
sequence, therefore, the two adjacent periods 57 and 58 are compared with the 
dynamic height currents. For station I the two intervals are 02J^0 to 07OO and 
0700 to 1120, and for station II O55O to 0750 and O75O to 1210 on 25 August. 
Of interest in this figure are the strong differences between the dynamic height 
current components and between the two sets of ^-hour averages that bracket the 
time interval when the dynamic height components were computed. 
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The results of these experiments, notably the results of l8 and 20 August, 
raise strong questions about the extrapolation of current data from any given 
single current meter station. There are often serious discrepancies between 
the ^-hour averages in corresponding periods between the two current meter sta- 
tions that are only one-half mile apart. Further, the drogues move near the 
current meter stations and show net transports over comparable ^-hour periods 
that are at variance with the average currents computed at the stations. 

Finally, the sequence of quarter-period current meter averages shown for 
the cases of the frontal passages on l8 and 21 August show that the dynamics of 
the response of currents to winds are not well understood. Moreover, in the 
period preceding the front of 21 August when the winds were from the southerly 
direction, the general tendency of the current meter data was to show an overall 
southward drift. 
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III. DETERMINATION OF EDDY VISCOSITY AND EDDY DIFFUSIVITY IN LAKE MICHIGAN 

Joseph C. Huang 

INTRODUCTION 



Theoretically, a flow field of incompressible Newtonian fluid can be 
solved by the equation of conservation of mass and by the Navier-Stokes equa- 
tion of motion, satisfying the respective boundary and/ or initial conditions. 
When the flow is viscous and laminar, a term appears in the equations of motion 
that is related to shear stress which, according to Stokes' law of viscosity, 
is 

du 
dy 

where the subscript "l" denotes ''laminar, " y is normal to the direction of the 
flow velocity, u, and ij, is the molecular viscosity (which is a property of the 
fluid depending only on temperature under ordinary conditions). In compressible 
laminar flow, when temperature variations are involved, the heat flux (following 
Fourier's law of thermal condition) is shown as 

^ dT 

where k is the molecular thermal conductivity, which is also a property of the 
fluid and usually depends only on temperature. 

On the other hand, if the flow is turbulent, a kind of motion with random 
irregular fluctuations is superimposed on the main flow. Since its fluctuations 
are irregular, it is impossible to describe the motion in all details as a func- 
tion of time and space. However, one can study its average behavior using the 
statistical methods pointed out by Hinze (1959). Therefore, in describing a 
turbulent flow in mathematical terms, it is convenient to separate it into a mean 
motion and a fluctuation, such as 



u == u + u 

T = T + T' 
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where the "bar" indicates mean quantities and the "prime" indicates the fluc- 
tuation from the mean. 

The presence of random irregular fluctuations of eddy motion of the flow 
appear as apparently increased viscosity or diffusivity (Defant I961) . In tur- 
bulent flow, therefore, additional Reynolds stresses must also be considered in 
the general governing equations of the motion. 

Analogous to the coefficients of molecular viscosity, Boussinesq (1877) 
introduced a mixing coefficient, A^, for the Reynold stress in turbulence ex- 
pressed in terms of the mean velocity gradients in the flow field, 

^ _ du 

t " " ^ u^v' " t dy 

where the subscript t denotes turbulent flow. In a similar definition to that 
for molecular diffusivity, Schlichting (1968) postulated that the turbulent 
heat flux is 

t P q qy 

The exchange coefficients of turbulent momentum and turbulent heat flux, 
A^ and A^, both have the dimension of a viscosity, \i (g/cm/sec). 

If the turbulence is isotropic (i.e., the relation between stresses and 
strains is independent of axis) then A. corresponds to the molecular viscosity 
\i, A corresponds to the molecular heat diffusivity (k/c ); and e = At/p, and 
the eddy kinematic viscosity corresponds to molecular kinematic viscosity 
V = [x/p* Turbulence shear stress is expressed as 

du 
T = p € — . 
t dy 

The eddy kinematic viscosity, e can be expected to occur as a constant only 
if the turbulence field is homogenous. Generally speaking, the eddy viscosity 
is 10^ to 10^ greater in magnitude than molecular viscosity, except in the vis- 
cous sublayer very near the solid boundary. The latter may be neglected in 
many cases to a good approximation. As pointed out by Lamb (1932), in describing 
the wind produced current, \i is replaced by A. . 

By the techniques mentioned above, it is possible to attack geophysical 
fluid problems which are usually turbulent in nature. In this way we can probably 
predict or evaluate the time-mean values of these complicated turbulences whose 
complete theoretical formulations have so far proved impossible (Lamb I952). 
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Most of the geofluid modelings are treated in this way. For example^ previous 
work by Sverdrup (19^7)^ Stommel (19^8)^ Munk (1950)^ Munk and Carrier (1950)^ 
Bryan (I965), and Birchfield (I967) all follow this kind of approach. 

Nevertheless not all the turbulence is homogeneous and isotropic. The 
eddy viscosity usually does have directional preference, Eddy viscosity be- 
comes effective only if there is some flow in the fluids and its value depends 
on the velocity distributions and on the characteristic length of the flow field. 
As shown in Boussinesq^s hypothesis and von Karman*s similarity theory 

t dy 'dy 

where L is the Prandtl's mixing length. Therefore A^ or e is not a property of 
the fluid itself^ whereas the molecular viscosity^ \i^ is. It is still extremely 
useful to introduce the eddy viscosity in relating the Reynolds stress to the 
mean velocity gradients of the flow. Though A, (or g) is not a physical con- 
stant^ but varies from one part of the flow to another, estimation of the magni- 
tude of the average eddy viscosity is still the key factor to the success of 
geophysical fluid modeling. In Lake Michigan, though the mean velocity is in 
general less than l/2 m/s (Ayers^ et al. 1958)^, the flows are turbulent in 
nature. Therefore, it is important to estimate the order of magnitude of the 
eddy viscosity and the eddy diffusivity in order to facilitate geophysical fluid 
modeling of Lake Michigan. 



THE MAGNITUDE OF EDDY VISCOSITY AND EDDY DIFFUSIVITY IN LAKE MICHIGAN 

In order to evaluate the magnitude of the eddy diffusivity and eddy vis- 
cosity associated with the thermal current structure of the late spring/early 
summer circulation in the Lake Michigan model, a series of experiments were 
carried out to verify order- of -magnitude values for the eddy diffusivity and 
eddy viscosity to be used with further theoretical investigations. 



An Estimate of Eddy Viscosity by Reynolds Stress and Vertical Velocity Gradients 

In August 1967^ a pair of current meter buoy stations were set 8 miles off- 
shore from Sheboygan, Wis. (see p. 12 of this report for a description of the 
buoy stations and their locations) » 

During the period of the current measurement, the research vessel INLAND 
SEAS was operating around the stations (during the daytime) whenever weather 
permitted. Wind, as well as other weather data, was automatically recorded on 
a digital meteorological tape recorder. The anemometer and vane that feed into 
the recorder are on the mast of the ISLAND SEAS ik m above the water surface. 
The Captain ^s log and hourly meteorological observations recorded wind from a 
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lower vane and anemometer that are about 11 m above the surface. Using the 
analyzed recorded mean wind versus height plotted on a log scale with the mean 
velocity zero at about 1 cm height (Floyd C. Elder, personal communication) the 
10-m level wind was interpolated with reference to the Captain *s log. The 10-m 
daily winds are shown in Table J-l. The vertical velocity gradients in the 
direction of surface current were deduced from two current meters separated 
vertically by 10 m. According to Deacon (I962), the Reynold stress is 



= A ^ 
t dy • 

— 3 

Then using p^ = 1 g/cm^ and p^^ = I.25 x 10 g/cm^ the eddy viscosity can be 
calculated. The data^ wind, wave conditions, absolute vertical velocity gra- 
dient, shear stress, and calculated eddy viscosity are all shown in Table 5-I. 
The values of Reynold stresses obtained here using Deacon's formula are in good 
agreement with the mean value of stress by Elder and Soo (1967) at The Univer- 
sity of Michigan Research Tower on the east side of the lake. 



An Estimate of the Upper Bound Value of Eddy Viscosity by the Result from 
Statistic Theory 

It is quite well known (e.g., Ayers, et al. I958) that a large clockwise 
eddy exists in southern Lake Michigan between Grand Haven and Michigan City. 
This large eddy has the length scale of a little less than half the dimension of 
the width of the lake. Bowden (196^1) has shown the effective eddy viscosity is 
proportional to the V5 power of the length scale, L, of the particular eddy 
under consideration. In Lake Michigan, let this largest eddy have average 
length scale of L = ^0 km or L = i+ x 10^ cm. Then, according to Richardson's 
formula, the eddy viscosity will be of the order lo'^ cm^/sec which is the upper 
bound of the numerical values of the eddy viscosity in Lake Michigan. 



An Estimation of Vertical Eddy Viscosity from Current Measurement Near the Bottom 

From 11 May to 11 October ±9^7 y a tripod- supported pendulum current meter 
was installed about 25O m offshore, T-l/2 miles south of Benton Harbor in a 
depth of 4.5 m. The pendulum was suspended only 20.3 cm above the sandy bottom. 
The monthly mean current is shown in Table 5-2. Using Lesser 's (195I) roughness 
length for sand bottom, z^ = .I5 cm, the friction velocity near the bottom can 
be determined through 

z + z 
\ = V[in(-^)]"' 
o 
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where u is the mean velocity, z is the height from the bottom during measure- 
ment, and Kq = O.Ul, The monthly mean friction velocity for z = 20.0 cm is 
shown in Table 5-2. The neutral stability condition generally persists near 
the bottom. Then the vertical eddy viscosity can be estimated as shown in 
Table 5-2 (Bowden I96J4) . 



TABLE 5-2. Estimated vertical eddy viscosity by current measure- 
ment near bottom. 



Vertical eddy 
Mean current Friction velocity viscosity 
(cm/sec) (cm/sec) (Az = k u z cm^/sec) 



Month 



o^ 



May- 


59.5 


June 


51.7 


July 


1+7.2 


Aug. 


58.5 


Sept. 


56.1 


Oct. 


69.8 



k.6k 

2.krj 
5.68 
I+.56 
U.58 
5.UU 



UO 

20.7 

50.9 
38.5 

56.0 
1+5.7 



Average value 



33 cnf /sec 



An Approximation of the Vertical Eddy Viscosity by Surface Wind 

Wind data from both the two- current-meter-buoy operation off Sheboygan 
and the wind and current station outside of the Benton Harbor beach are used 
for approximating the vertical eddy viscosity in Lake Michigan. In the two- 
current-meter-buoy operation off Sheboygan, wind data were measured aboard the 
R/V inland seas. At the current meter station outside of the Benton Harbor 
beach, the wind was measured Just above the water. From Sverdrup, et al. (I96I+), 
depending on whether the magnitude of wind was greater or less than 6 m/s, the 
vertical eddy viscosity was approximated as shown in Table 3-3. Notice that 
the vertical eddy viscosity is smaller when the wind is averaged for a longer 
period of time . 
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TABLE 3-3. Estimated vertical viscosity by wind speed. 

Vertical eddy 
Time Wind (m/s) viscosity, A^ cm^/sec Remarks 



IT Aug. 


1967 


9.5 


18 Aug. 


1967 


9.8 


20 Aug. 


1967 


5.7 


21 Aug. 


1967 


5.^ 


25 Aug. 


1967 


6.2 


2k Aug. 


1967 


h.6 



581 

165 
100 



Wind averaged 
over 6 hours 



Average value 



206 cm^/sec 



May 1967 


3.5 


June 1967 


2.5 


July 1967 


5.1 


Aug. 1967 


5.2 


Sept. 1967 


2.8 


Oct. 1967 


k.l 


Nov. 1967 


k.k 



i;^ 


Monthly mean 




wind 


16 




50 




55 




22 




TO 




87 





Average value 



^3 cm^/sec 



A Measurement of the Eddy Diffusivity by the Dispersion of Marked Particles 

In a series of experiments conducted in Lake Michigan, the fluorescent dye 
rhodamine B was used as a tracer for diffusivity measurements. Rhodamine B has 
peak ultraviolet fluorescence wavelength of 58O millimicrons. In May I968, a 
fluorometer with only one sample cell was used. Because in this season the 
phytoplankton was so abundant, and since they have an ultraviolet fluorescence 
frequency Just a little lower than that of rhodamine B, the fluorescence of the 
peak concentration of rhodamine B dye was not distinguishable from the background 
after an hour or so. In all subsequent experiments, a comparative measuring 
fluorometer with two cells was used. The rhodamine B mixed with natural lake 
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water flows through the dye sample cell, while the standard cell contains natural 
lake water. The relative difference in fluorescence between the two cells is 
measured in order to obtain reliable data in the possible presence of any inten- 
sity variation of ultraviolet light lamp,^ systematic changes in the characteris- 
tics of photocell, and differences in the amount of naturally occurring fluorescent 
substances in the water (Noble and Ayers I961) • Both instantaneous source and 
continuous source measurements have been made in I968. 

Instantaneous Sources . Three runs of instantaneous sources of marked 
particles for diffusivity studies were conducted. These were 1 and 8 miles off 
Grand Haven in May, and 5 miles off Milwaukee in June. Rhodamine B dye was 
diluted with methanol to have a specific gravity very close to the lake water. 
Generally, one part of stock solution of rhodamine B is Uo^ acetic acid solu- 
tion was mixed with 5-1/^ parts of methanol. The dye mixture was gently poured 
out on the surface lake water at the upwind side of the boat in order to allow 
the boat to drift away with least interference to the dye patch. The first and 
second runs were investigated aboard the r/v MYSIS and the third run aboard the 
R/V INIAKD seas. During the first 20 minutes after the source had been instan- 
taneously releaseid, the diameters of the expanding dye patch were estimated by 
eye with the ship's length as reference. At 20-minute intervals, the boat was 
used to coast through the dye patch longitudinally and laterally at the lowest 
possible speed. For each transect, fixes were obtained from shore objects when 
entering and leaving the patch, along with the time required to traverse through 
the patch. Continuously running back and forth across the dye patch, the con- 
centration of the dye patch was continuously recorded from the fluorometer on an 
AZAR recorder. The concentration distribution of a typical pass through the 
patch is shown in Figure J"!- Most generally, the dye concentration distribution 
showed a form of Gaussian distribution. The diameter of the patch, longitudinally 
or laterally, is properly approximated as k times the standard deviation of the 
patch from its mean position. Batchelor (1955) shows that (assuming the fluid 
is in a uniform, homogeneous, and steady, mean velocity field) at the initial 
stage of diffusion process, the effective diffusivity (which is proportional to 
the rate of change of the square of standard deviation with time) is proportional 
to the first power of the elapsed time; in the intermediate phase, it grows as a 
V5 power of the patch size; and at the final stage, as the standard deviation 
grows as l/2 power of the elapsed time, the diffusivity should asymptotically 
approach a constant. The general information and the value of effective diffu- 
sivity are shown in Table J-U. The growth of dye patch versus time is shown in 
Figures 5-2, and 5-5. 

Continuous Source . Four runs of continuous (in time) point sources of marked 
particles were conducted near the meteorological research tower 1 mile off the 
east shore in Lake Michigan. Ten gallons of rhodamine B and methanol mixture 
were contained in two 5-gallon carboys connected together by plastic tubings as 
shown in Figure 5-^. The continuous point source was released from the apparatus 
through a rubber tubing which had its outlet about 20 cm above the water surface. 
The rate of discharge of the dye was controlled by a screw clamp on the tubing,, and 
was usually allowed to flow at a rate of 5 gallons/hour, or approximately h g/sec. 
In a homogeneous fluid with a steady uniform current, the continuous dye source 
was observed as a slender plume extending downstream. 
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Figure 5-1- A typical distribution of dye concentration across a dye patch. 



TABLE 5-i^■. General information and effective eddy diffusivity ©f dye 
patch studies. 



Run 


Run 1 


Run 2 


Run 5 


Date 


1 May 1968 


2 May I968 


27 June 1968 


Location 


1 mile offshore 


8 miles offshore 


k miles offshore 




near Grand Haven 


near Grand Haven 


near Milwaukee 


Depth 


20 ft 


ko ft 


65 ft 


Duration of run 


1150 - 1?10 EST 


1500 - IU50 EST 


0910 - 12if0 EST 


Current velocity 








cm/sec 


11.6 


5.5 


7.6 


Wind 








Direction 


m 


SW 


SSW 


Speed (kts) 


15 


15-20 


5 


Lake condition 


2-5 ft waves 


1 ft wave 


2 ft waves 


Effective eddy dif- 








fusivity (cm^/sec) 








Lateral 


Not available 


lli^O 


7^5 


Longitudinal 


Kot available 


^131 


^-019 


Remark 


Very poor visi- 
bility 
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Figure 3-2. Lateral growth of the dye patch with time, 
(a) 2 May I968. (b) 27 June 1968 
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Air Intake- 




Screw Clapper* 



Figiire 5-^. Continuous source apparatus. 



The direction of the dye plume was reported by the observer using a theodo- 
lite mounted on the meteorological research tower 57 ft above the water surface. 
One or two buoys were dropped as markers on the edges of the dye plume at right 
angles to the direction of flow. The appropriate distance of the marker from 
the source depended on the mean velocity of the current which was estimated from 
the drogue study set earlier. 

The r/v MYSIS crossed the dye plume back and forth at 10-minute intervals 
at the lowest speed she could make (about k.^ mile/hour). As soon as the dye 
plume passed through the markers, samples of water were taken through a constant- 
flow fluorometer that fed into an AZAE recorder. Communication between the 
tower and the boat was maintained through a VHF transceiver. Whenever the boat 
entered and left the dye plume a "MARK" signal was transmitted to the theodolite 
observer. The observer obtained the azimuths of both edges of the plume and the 
elevation from the tower to the boat. 

The horizontal distance between the boat and the tower was obtained from 
the elevation angle as well as by fixes from land objects. The angle subtended 
by the edges of the plume is the difference of the two azimuths. The arc length 
of the plume at this distance from the tower was arbitrarily taken as four times 
the lateral standard deviation unit. At each return crossing, the vertical 
diffusion was measured by lowering the fluorometer intake at half -meter inter- 
vals until there was no indication of the dye. Half of the maximum vertical 
displacement was taken as the vertical standard deviation unit. 
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At the last few minutes before the exhaust of the point source^ continuous 
crossings were conducted at various distances along the plume from the tower 
to approximately I-I/2 miles downstream. These crossings measured the asympto- 
tically steady diffusion of the plume. 

These continuous point source diffusion studies were carried out on 50 July 
and 1-2 August I967. On 50 July, sampling was discontinued about 2 hours after 
the experiment was initiated, because of high seas. The dye plume was contin- 
uously sampled for 5 hours on 1 August, and during the sampling period the cur- 
rent shifted its direction of flow about I5 degrees. Two sampling runs were 
conducted on 2 August under a very steady current condition. The general infor- 
mation for each run and the value of diffusivity are shown in Table 5-5. The 
lateral growth of dye plume versus time at a constant distance from source is 
shown in Figure 5-5. The lateral growth versus distance under asymptotically 
steady condition is shown in Figure 5-6. 



A Measurement of the Horizontal Diffusion by Drogue Studies 

A drogue is made of k pieces of light sheet metal (e.g., aluminum alloy) 
approximately 1 x 2 m in dimension, linked together in the shape of a cross 
and suspended by chains in the water by a float which supports a radar antenna 
above the surface (similar to Farlow I965). It is frequently used for measuring 
the mean current of the flowing layer at the drogue depth. The equivalent radius 
of the drogue is about 1 m and, based on its geometry and the corresponding 
Reynold numbers (10^ to 10^ in the Great Lakes), its time constant is estimated 
to be 100 sec.(Okubo and Farlow I967) . Though the small-scale turbulence may 
be averaged out, the drogue is adequate for measurement of large-scale turbulent 
diffusion for eddies larger than 10 m in scale which represent the energy- 
containing eddy scales in the Great Lakes (Csanady I965; Okubo and Farlow 1967) . 

On 27 June I968 a hexagonal pattern of surface drogues, with the center 
one very close to a fixed buoy, was set about 5 miles offshore near Milwaukee 
in Lake Michigan. The original distance between each pair of drogues was l/2 
mile. After setting the pattern, the r/v INLAND SEAS continuously fixed each 
drogue for the following 6-I/2 hours. Similar to the result mentioned by 
Okubo and Farlow (1967), the growth of the pattern size with time was not 
clearly noticed. The diffusion constant was deduced from the second moments 
about the mean of the haxagonal pattern. At any particular time, the instan- 
taneous positions of all drogues were interpolated from their respective pre- 
vious fixes. The general information on the experimental run and the value of 
eddy diffusivity are shown in Table 5-6. 



COMPARISON OF THE VALUES OF EDDY VISCOSITY AND EDDY DIFFUSIVITY WITH OTHER 
STUDIES 

In Lake Michigan, there is little data about eddy viscosity or eddy dif- 
fusivity. Some figures from the ocean or other lakes are available as reference « 
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Figure 5-5. Lateral growth of dye plume at a constant distance 
from the source, (a) 1 August I968. (b) 2 August I968. 
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TABLE 5-6. General information and eddy diffusivity of 
drogues study in Lake Michigan. 



Date: 

Location: 

Depth of water: 

Depth of drogues: 

Duration of run: 

Current velocity: 

Wind: 

Direction 
Speed 

Lake condition: 

Mean initial separation: 

Final deviation from mean: 

Effective eddy viscosity: 



27 June 1968 

5 miles off Milwaukee 

65 ft 

5 ft 

I5OT - 1850 EST 

7.6 cm/sec 



N 

6 kts 






1/2- 


1 ft 


waves 


2578 


ft 




3UOO 


ft 




5A X 


: 10* 


cm^/sec 



In modeling the Gulf Stream, Munk (1950) estimated the lateral eddy viscosity 
of 5 X 10'^ cm^/sec. Stommel (1955) estimated the eddy viscosity of Gulf Stream 
near the Florida Straits to be 2 x 10^ cm^/sec. Defant (I96I) determined the 
exchange coefficient of lateral large-scale turbulence in ocean phenomena as 
5 x 10'^ cm^/sec. 

In general, the vertical turbulence scale is much smaller in dimensions 
than that of the horizontal. Therefore the magnitude of vertical eddy viscosity 
is much smaller. According to Defant (1961), Sverdrup obtained the value of 
the vertical eddy viscosity of 75-260 cm^/sec in the North Siberian Shelf -re- 
gion; Fjeldstad estimated the vertical eddy viscosity in the Caspian Sea of 
0-22if cm^/sec, and Suda measured the value of 680-7500 cm^/sec for the Kuroshio 
and of 150-1^^60 cm^/sec in the Japan Sea. 

For vertical eddy diffusivity, some typical values have been summarized by 
Defant (I961) as: Mediterranean, k2 cm^/sec; California current, 50-^0 cm^/sec; 
Caspian Sea, 2-l6 cm^/sec; Equatorial Atlantic Ocean, 520 cm^/sec. Csanady 
(1965,196^,1966) has measured the eddy diffusivity in Lake Huron and Lake Erie, 
and found horizontal diffusivity of the order of 5OO-25OO cm^/sec in Lake Huron 



k2 



and, 1000-5000 cm^/sec in Lake Erie. These values are much greater than verti- 
cal eddy diffusivity, which is only of the order 0.1-10 cm^/sec. Okubo and 
Farlow (1967) measured the horizontal diffusivity for large eddies in Lake 
Michigan and Lake Erie using drogues and found the eddy diffusivity of 5 x 10 
to 6 X lO''^ cm^/sec. 

According to our estimation, the vertical eddy viscosity is on the order of 
55 to 2 X 10^ cm^/sec in Lake Michigan. If the vertical eddy viscosity is 
usually two orders smaller than the horizontal value (Defant 196l)jj then the 
horizontal eddy viscosity should be in the range of 3.5 x 10^ to 2 x 10^ cm^/sec, 
which is within the upper limit of lO"^ cm^/sec as we have estimated. The 
lateral eddy diffusivity we found using marked particles is 750-1200 cm^/sec for 
small eddies and ^.k x 10^ cm^/sec for larger eddies. The effective longitudinal 
eddy diffusivity for small eddies is in the order of ^ x 10^ cm^/sec. All the 
data we used appear quite consistent and in good agreement with that of others. 



CONCLUSION 

In modelling a geofluid problem^ the eddy viscosity and the eddy diffu- 
sivity are of critical importance in order to predict the natural current or 
wave phenomena with a similarity to laminar flow. In Lake Michigan, though the 
mean current velocity is small in general, the flow field is turbulent in 
nature. In solving or explaining the flow pattern of the mean lake current, it 
is possible to use the governing equations of laminar flow with eddy viscosity 
and eddy diffusivity in place of the molecular viscosity and the molecular dif- 
fusivity* A series of experiments have been conducted in order to estimate 
these two values. The data presented are quite consistent and in good agree- 
ment with the data reported by other investigators. The vertical eddy viscosity 
in Lake Michigan is in the range of 1 to 10^ with a mean 10 cm^/sec. The 
horizontal viscosity is in the range of 10^ to 10"^ with a mean value of 10^ 
cm^/sec. The eddy diffusivity may reach the same magnitude as the viscosity 
but it is in general smaller. A typical mean value for vertical eddy diffu- 
sivity is 5 cm^/sec^ and for the horizontal eddy diffusivity 10^ cm^/sec. 
Though few experiments are far from enough to state the characteristic values 
of turbulent phenomena, the data presented together with values reported by 
others can be used as good reference values for Lake Michigan. 
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IV. MJMERICAL EVALUATION OF STERN ^S CIRCULATION MODEL 

James H. Saylor 



DEVELOPMENT OF STERN^S MODEL 

In 1965^ Stern presented a model for the interaction of a uniform wind 
stress with a geostrophic vortex. The results of these calculations showed 
that a geostrophic vortex subjected to a uniform wind stress would move, as an 
entity, with an Ekman drift, provided that the thermocline was "soft" and that 
internal waves could be generated at the center of the vortex. Stern* s results 
have been evaluated in terms of wind stresses and vortex dimensions that might 
be expected to occur in Lake Michigan (Noble 196?). The results predict a very 
broad range of internal wave periods that are within the limits of values 
measured for Lake Michigan (Mortimer 1965). 

The starting point for the theoretical model is the assumption of a hydro- 
static, two-layer model; upper layer of density p and lower layer of density p^; 
no friction; and the fluid incompressible. 
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Dv J ^ VP 


A J1 


— - + fkxv+ — = • 


-B'^*^ 


Dt p 



(1) 



where H = depth of upper layer. 



and 



p - p 
o 



D d ^ 

^ A A A A A 

U = V + Wk = Ui + vj + Wk 



A 

V-U = 



Let f/h = (S/Sz) G where 9 is permitted to be a lateral wind stress only. ^In par- 

^ / \ A ^ 

ticular, take (x,y,z = o,t) equal to the surface wind stress, t. Let U = 
Uq + U-j^, where 11^ is taken as the geostrophic component of (l). 
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A ^ h 

U = V + KW 
o o o 



A / A 

u = V + m^ 

b b b 



Thus 
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A /I VP / 

+ fk X V + — = - g'k (2) 
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where 



A 

V-U =0 , 
o 

A A 

and whence^ U-j^ is readily defined as the component of U due to the surface wind 
stress, T. Call this part the "Ekman drift." Boundary conditions are established 
by requiring no displacement of the free surface^ i,e,^ ^'U^ = - k^^U^ = 
W^(x^y^z = O)^ and by requiring at the fluid interface that the "Ekman drift" is 
limited to the upper layers^ and that the wind current joins the geostrophic 
component smoothly Just above the interface, 

^A 

U^(x,y,z = -H,t) =-^(x,y,z = -H,t) = , 



which implies 



5. ..^ ^ ^^ n 

9(x,y,z = -H,t) = ^ = ^ = 



A/ s A 

9(x,y,z = 0,t) = T 



Subtracting the geostrophic component, equation (2), from the momentum equation, 
equation (l), gives: 

A 

which can be written as: 

b^AAAAA Aa^G ,, 

-T— + V • Vv + V -Vv + V -Vv + f k x V = -:r- . ( 5 M 

at boobbb bdz ^^ ^ 
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Taking the vertical component of the curl of equation (5')^ denoting: 

if s k.V X U = k-V X V 
^o o o 

A 4 4 A 



Stern obtains, 



where: 



^ A /I /I A 

2^ = k-V X U = k*V X V 

^b b b 



(f + <f )-^— + -^(k-V X 0) = v^-ViT + Z T. (4) 

o' az 3z' b ^o i ^ ^ 

1=1 









dz dx o b dz dry' o b 

o D o b 
Ts 



Sz 5x dz ^ ' 



a 



<^ 



^ .. ^f 



Ts = (W + W^)^— + W^ . 
o b dz b dz 

Integrating equation (h) between z = -H and the free surface: 

/^ (f + ^ )^ dz + /^ |-(k-V X 0)dz = /^ v^.V^'dz + /^ Z T.dz (5) 

-H ^ ^o^Sz -H Sz^ ^ -H b o -H . , i ^^^ 

1=1 



which becomes: 



•(f + W )W^ - k-V X T + f V -V? dz = - / Z T.dz . (6) 

o -H b o -H . , 1 

1=1 



i+8 



stern compares the magnitude of the terms in equation (6) and shows that if the 
Rossby numbers based on both the geostrophic and "Ekman" currents are small, 
i.e. y 

D ^ O _ 

the / T. terms are of negligible magnitude. Defining the wind-driven volume 
transport as; 



/^ v^dz = M = - - k X T 
(which is the Ekman transport in the conventional sense). 

(f + ^ )W^ - k-V X T + M-V^ = (6') 

where A a h 

^ ^^ k X T-V6 

To include the effect of f varying with latitude, use the beta-plane approxi- 
mation, which simply replaces the V<^ term with V(f + ^ ). So, the general so- 
lution is written: 



which may be expressed as: 



«* = -'-|ff|^ (T") 



where 



For the vorticity of the geostrophic component in the upper layer we have the 
relation: 



^9 






Integrating this expression between -H and the surface: 

,0d,. .,, ,0 ^ 



Ih H(' ^ i'- = /^ S^ - ' 






A A 

dt 



which may be written as: 

d(f + ^ ) 



dt H 



, f rr ^ X T f dH 



f + £' ~ H dt 



(9) 



Equation (9) is the generalized equation of motion for the geostrophic component 
in the mixed "Ekman layer." 

For a constant, uniform wind stress (v x r = 0) and for a limited area in 
which f ^ constant, equation (?') reduces to: 

(k X T)*Vf 
^^=-77— F^ (10) 



^^-o 



o 



and the kinematic equation, equation (9)i becomes: 

H dt " H ( f + f' )^ 



dt 



. (11) 



Stern's analysis with the necessary conditions and restrictions, describes the 
interaction of a localized, geostrophically-balanced flow (i.e., within the 
confines of the simple, two-layer model envisaged) with a wind-driven "Ekman" 
current. The vertical velocity (W^) generated by this interaction is propor- 
tional to the gradient of the geostrophic vorticity (V^), a quantity which is 
conserved in the upper layer, and to the wind stress, t. In this manner, it 
is seen how the vertical velocity transfers momentum across the interface to 
the lower layer, and how there is a tendency for the vorticity of the water 
column to increase due to the advection of neighboring vorticity by the mean 
Ekman drift. 
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APPLICATION OF STERN'S MODEL TO LAKE MICHIGAN 

If the wind stress (t) and the coriolis parameter (f) are independent of 
position, the last term of equation (ll) reduces to: 

^(kXT)-V^^ ^k.(V^^XT) k.Vx(^^fT) 

a k.v X (— ) 

which is equivalent to a force per unit mass of - ^^r/Ef acting on the upper 
layer. 

Consider a two-layer, model (Fig. k-l) with Ekman drift at right angles to 
T and directed along the x - axis, i.e.. 



T = positive constant 

y 

T = 0. 

X 

^ ^_^ ^ ^ ^ ^ -^ ^ ^ ^ - • ^ ^ X^ 

Vorticity <, I E^ /> .y -\ + \ 



/// ^ /////// //> ^ ' 



Vorticity li a. ^ H^ '• ^ ' 5- -H, - H, 

Figure i|.-l. Schematic two- layer model. 



Let the positive displacement of the interface at z = -Hi be small and given 
"by Ti(x,y, t). The vorticity equation for the lower layer is: 

d4 _ ^ ^ 

dt az • 



Integrating between z = -H^ - Hi and z = -Hi + r| 



,-Hi+Ti d^ _ pHi+n ^ ^ ^^ 

iHi-H2 dt ^^ iHi-He ^ hz ^^ ' 
dt Ha' 'z = -Hi+Ti , 
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^ ^z = -Hi + n dt 



Integration gives: 



^*#- = ^*it" ' 



t=0 t=0 



^^(t) = Ul^llll . |^^(^^y^t=0) . (12) 



The vorticity equation for the upper layer is: 



'i^^i'^^^^t'" 



where U = r/fHi = mean "Ekman" drift. 



(15) 



The interaction of an Ekman drift current with a geostrophic vortex can be 
described in terms of a vertical suction velocity W^, which is proportional to 
the gradient of the geostrophic vorticity and to the wind stress. This interac- 
tion leads to an increase in the total vorticity of the water column with time 
because of the transfer of momentum across the interface. It was shown that if 
T and f were constant, the result of the interaction can be described as advec- 
tion of vorticity at the rate of mean Ekman drift. Equation (I5) expresses this 
fact. If the geostrophic vortex is centered at x = y = at time t = 0, at 
which time an Ekman current is established, equation (I5) states that part of 
the vorticity will stay at x = y = (represented by the interface disturbance), 
while part is advected with the drift current. 

An expression for the vorticity difference between the upper and lower 
layers comes from Mar gules * formula 

*ri =<& - f^vit) (Ik) 



where 



r^ -^ ^ 

^ -^""W 



From (12), (15) and (lU): 



/l. + ui-)(^ .^Ivin) + — -^ = U — ^ 



^ ' "^ hi'^-YLz f '^"'^ ' Hi ^ "" H2 ax ' ^"""^^ 
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Separating equation (l4) into steady-state and time- dependent components by- 
setting 



Tl(x,y,t) = Ti^(x^y) + 7i^(x,y,t) , 



for the steady- state component: 






for the time- dependent component: 

Assuming a solution of the form: 

ik(x-Ct)+i% 
''t " ° 

(C = phase speed) then 



\ = ^ 



vin^ = -(k- + i^)!^^ 



equation (I7) becomes: 

^^^Hi_^^ + g' p(k^ + -«^)] = U[l H- g' ^(k^ + i^)] , (17.) 



or, 



§=[l+|(l-k--^0"" (18) 



where 



K = (^)'/^ k 



^, = (^)V^ i 
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-4 -1 

Take f = 10 sec (an approximate value for Lake Michigan). 

If k^ » Hg/Hi, C 21 U and short waves are advected at the speed of the 
mean Ekman drift. 



If k^ = i^ is small- the phase speed is: 

(1 + k^ + 4)Hi 

and a surface vortex of radius ^(g* -3)^/^ will leave the origin in time: 

^"2^ f^ ^ ' LuHi(l + 4 + 4) J • 

To estimate the eddy diameters, ratio of propagation, and internal wave 
periods involved, use the following parameters as indicated in Figure i|— 2. 



I - ° 



a 



j^ = 20 m ^ = 0.9982 (20° C) 






H2 - 100 m /tf© s^+A^= 0.9998 (8" C) 

Figure k-2. Specific parameters of the two- layer model. 



Let: 



^:^ = ^:^ ^ f = 10" sec" 



:' = (— )g « 1.57 X 10"^ m/sec^ 



T = 2 X 10' (m/sec)^ 



Then, 



*T* 

U = — - = 10 cm/sec 
rHi 



% 



^ ~ ^[{1 H- k^ + 4)Hi + hJ 
eddy diameter = — ( ^g^ )"'"/^ 



TABLE i+-l. Calculated:eddy diameters, phase velocities, 
and interval wave periods of upper layer 
in Lake Michigan. 



K = h 


Eddy 

diameter, 

km 


cm/sec 


Phase speed 

c, 

cm/ sec 


Wave 

period, 

days 


1 


59 


10 


5.75 


12.1 


2 


1.95 


10 


5.0 


^.57 


5 


15 


10 


6.67 


2.25 


it- 


7.8 


10 


9.1 


0.99 


5 


^.8t 


10 


9.55 


0.59 



Effect of Density Difference, Ap 

Keeping other parameters fixed, let Ap = l/2 of former value, i.e., take 
g' = 0.8 X 10 (cm/sec^). Considering the density change, the result is cal- 
culated as shown in Table 4-2. 



TABLE 4-2. Calculated eddy diameters, phase velo- 
cities, and internal wave periods for 
upper layer in Lake Michigan. 



Eddy 
k = i diameter, U, 
km cm/ sec 



Phase speed 
cm/sec 



Wave 

period, 

days 



1 


28.1 


10 


3.75 


8.7 


5 


9A 


10 


6.67 


1.63 


8 


5.5 


10 


9.55 


O.J+5 
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Effect of T 



The wind stress determines the speed of the Ekman drift. Within the limits 
of (v-[^|fL « 1 then^ the only effect of t is to determine the rate of advection* 
If the parameters of the basic model remain fixed, and t is reduced to l/2 of 
its former value (1 x 10 (m/sec)^), the solutions are modified as shown in 
Table k-3. 



TABLE ^4-5. Calculated eddy diameters, phase velo- 
cities, and internal wave periods for 
upper layer with wind stress changed. 



Eddy 



^^ ~ ^^ diameter, 
km 



Phase speed Wave 
U, C, period, 
cm/sec cm/sec days 



1 


59 


5 


1.87 


24.2 


5 


15 


5 


5.55 


i^.5 


8 


k.Qj 


5 


U.76 


1.18 



Effect of Depth of Thermocline 

Assume the same parameters and overall depth as basic model, but let Hi 
10 m, H2 = 110 m, then the results will be as shown in Table h-k. 



TABLE k-k. Calculated eddy diameters, phase velo- 
cities, and internal wave periods in Lake 
Michigan with the depth of thermocline 
changed. 



K = h 


Eddy 

diameter, 

km 


cm/sec 


Phase speed 

c, 

cm/sec 


Wave 

period, 

days 


1 
5 


i^5.4 
14. U 


20 
20 1 


4.28 
12.6 


11. T 
1.52 



For Hi + H2 = constant, the effect of decreasing Hi is to increase the 
mean Ekman transport velocity (u), to increase the wavelength associated with 
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each wave number^ and to decrease the wave period (because of increased advec- 
tion rate) . 

The range of values for internal wave periods calculated from this model 
are within those reported by Mortimer (I965) and also within the range of perio- 
dicities of current speed fluctuations calculated by Verber (1965). However, 
because of the limitations imposed by the assumptions of a simplified model, 
these values cannot be considered to be definitive. 

Stern* s model does provide mechanisms for the generation of phenomena 
observed in field measurements, such as rotary currents from buoy stations, 
rotary drogue motions, internal waves, and a transfer of wind energy into the 
water circulation in a way that might explain the buoy station measurements re- 
ported in the previous section of this report. 

It will be necessary to design future field experiments relating the three- 
dimensional temperature field to the three-dimensional current field as a func- 
tion of time. These experiments will be time-constiming and expensive, but the 
numerical results presented in this paper indicate that there is a strong reason 
to make definitive tests of the Stern model. 
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V- THERMAL CURRENT STRUCTURE IN LAKE MICHIGAIT, ASSOCIATED WITH 
THE SPRING WARM-UP SEASON— A THEORETICAL STUDY 

Joseph C. Huang 



INTRODUCTION 

The early spring condition of the lake lends itself most readily to mathe- 
matical analysis. During the winter^ the water mass of the lake becomes nearly 
isothermal^ and has a three-dimensional circulation pattern reaching to the bot- 
tom of the lake basin. In May^ "while the main water mass is still isothermal 
and at its winter temperature of ^'^C or less^ the edge of the lake begins to 
warm under the influence of a combination of factors, including warm-water run- 
off from the adjacent land. The nearshore boundary water warms to a temperature 
about S^'C (and is isothermal to the bottom) in water depth of 20 m, and extends 
a mile or more offshore. The edge between the boundary water and the main water 
mass of the lake is characterized by a sharp temperature gradient (the "thermal 
bar"), water color differences, marked differences in phytoplankton populations, 
and strong geostrophic currents along the temperature gradient. As the season 
progresses, the boundary water edge becomes more complex and often shows multiple 
bands of water color and temperature 2 to ^ miles or even more offshore, parallel- 
ing the perimeter of the basin. The period of 28 May to 5 June 1967^ when this 
process operated, has been described for Lake Michigan by Stoermer (I968) and 
by Noble and Anderson (1968). This same process has been well documented for 
Lake Ontario and Lake Huron by Rodgers (1965,1966). 

In the Great Lakes area, the wind system is very much modified by the con- 
tinuous presence of a high pressure in the spring- summer period and a low pres- 
sure in the fall -winter period due to the meteorological effect of this world ^s 
greatest freshwater masso The wind observed in the lake is quite different from 
that of nearby land stations o Both the magnitude and the direction of wind are 
frequently changing, seldom consistent over a long period of time^ The time- 
averaged net energy put into the lake by wind stress may not appear as a simple 
wind-current relationship. This is especially obvious during the lake warmup 
season. In the spring, when the lake water is still relatively much colder than 
the already warmed adjacent land, the air mass above the lake is frequently in 
a stable condition. Wiresonde cross- sections in Lake Michigan have demonstrated 
the intense springtime inversion prevalent before the warming of the water sur- 
face (Bellaire I965) . At the Wisconsin shore, it has been observed that water 
remains smooth even beneath offshore winds of greater than 10 knots (at 50 ft). 
However, as far as the thermal energy connected with temperature difference be- 
tween land and lake and air and water are concerned, temperature plays a much 
more constant role than the varying wind all the year around. During the lake 
warmup season, temperature gradients, either horizontal or vertical, show per- 
sistence, which, in turn, suggests that the thermal energy input may dominate 
the lake circulation dynamics, 
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Noble (1967) has indicated that the circulation dynamics of Lake Michigan 
are probably dominated by geostrophic forces rather than the wind stress during 
the summer season. It is the purpose of the present study to model Lake Michi- 
gan, considering the geostrophic thermal gradient force as the prime source for 
the lake current generation, and to investigate the thermal circulation pattern 
of the lake as a comparison with the general current reported by others (e.g., 
Ayers et al. 1958). The study has not reached the final stage yet. The model 
may be refined by some modification much closer to the real simulation of the 
lake. The preliminary result shows that the thermal circulation pattern in 
Lake Michigan is in very good agreement with general current data reported 
earlier (Ayers et al. 1958; Rodgers 1965) as well as our own investigations. 



MODELING CONSIDERATIONS AND ASSUMPTIONS 

In the geophysical fluid model of large scale motion, Coriolis effect must 
be taken into consideration. During the spring warming season, a sharp hori- 
zontal temperature gradient always exists within a few miles from shore around 
the perimeter of the lake. In a stable atmosphere, when winds are constantly 
shifting above the lake, the thermal gradient force comes to equilibrium with 
the Coriolis force. Generally speaking, only the geostrophic thermal gradient 
system is apparently in operation during the spring warmup season. Equilibrium 
between the temperature gradient force and the Coriolis acceleration produces a 
component of horizontal velocity perpendicular to the direction of the imposed 
thermal gradient in a rotating system. This is the geostrophic thermal gradient 
relationship which is considered dominating. 

In Lake Michigan, the length (about 500 km) is much greater than the width 
(mean about 120 km). The wide Straits of Mackinac connect the northern end of 
Lake Michigan to an even larger lake — Lake Huron. In the Lake Michigan model, 
to the lower order of approximation, the flow field, away from the north and 
the south ends, can be considered as zonally independent. 

The Coriolis force varies with geographical latitude. In Lake Michigan, 
its magnitude ranges from 0.^88 x 10 rad/sec in the southern end to 0.525 x 
10 rad/sec in the northern end. The variation from south end to north end is 
small. The Beta- effect in the Lake Michigan model can be neglected. 

An east-west cross section of the central part of Lake Michigan is shown 
in Figure 5-1. The deepest point in Lake Michigan is 28I m in the northern 
half of the lake. The southern portion has a maximum depth of about I50 m. 
In the Lake Michigan model, with a width of 120 km and depth of I50 m, a flat- 
bottomed trapezoidal cross-section as shown in Figure 5-2 is not too far from 
reality. Since the dimensional reference length, the width of the lake, is 
small compared to the radius of the earth, the centrifugal acceleration can be 
neglected. 
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Figure 5-2. Cross-section of a long, symmetrical trapezoidal 
Lake Michigan model. 



During the spring warming season, though the temperature structure in Lake 
Michigan is in a continual changing semi- steady state, the rate of heating in- 
creases with only some slight change in magnitude. Therefore the thermal body 
force is a seasonal function varying slowly with time. In modelling the lake 
with respect to the warmup season, all variables can be considered as time 
independent without losing the good approximation. 

The frictional forces are important near all boundaries. In the flow field 
of a turbulent nature as in Lake Michigan, the frictional forces are expressed 
as the space rate of change of the products of eddy viscosity and Reynolds 
stresses in analogy to the laminar flow. The average horizontal eddy viscosity 
in Lake Michigan is the order of 10^-10^ cm^/sec and the vertical eddy vis- 
cosity (1-10^ cm^/sec) is much smaller. Eddy diffusivity may be of the same 
order as eddy viscosity but it is a bit smaller in general. In our preliminary 
model 10^ cm^/sec is used as the average value for eddy viscosity and eddy 
diffusivity in both horizontal and vertical direction for ease in management. 
In the future refinement, separate values for horizontal and vertical viscosity 
are necessary in the Lake Michigan model. Eddy viscosity and eddy diffusivity 
of the fluid depend not only on temperature (such as the molecular viscosity), 
but change from place to place and time to time, depending on the local charac- 
ter and dimension of the flow. The overall average values of the order of mag- 
nitude of eddy viscosity and eddy diffusivity can be treated as "constant" for 
application to the geophysical fluid problem (see section III of this report). 

In general, though a weak thermocline exists during the later portion of 
the warming-up season, the lake is homogeneous in this particular period. In 
our preliminary model, we shall treat the fluid as homogeneous. 

Boussinesq approximations are used in the formulation, and the equation of 
state is a linear relation between the density and temperature as 



p' - P^[l - a(T' - T^)], 
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where a. and T are specified in Table 5-1. Furthermore, the dissipation func- 
tion is neglected in the heat equation. 



TABLE 5-l» Nomenclatures and nondimensional parameters. 



Symbol 



Definition 



X 
Y 
Z 

q. 

p 



V 



AAA 



Ot 

A 
Aq 

a 

L 

D 

AT 



Horizontal east-west coordinate, meridional direction 

Symmetrical north- south coordinate, zonal direction 

Vertical coordinate, negative is in the direction of gravity 

Velocity with components (u^v,w) 

Deviation of the pressure from hydrostatic pressure 

Temperature 

Constant lowest temperature near the bottom at the center of the 
lake 

Two-dimensional gradient operator, i -r- + k -r- 

ox oz 



(N.B. The above symbols primed have dimensions, unprimed are dimen- 



Two-dimensional Laplace operator, 
(N.B. The 
sionless. ) 

Unit vector in the (X,Y,Z) directions 

Acceleration due to gravity 

Coefficient of thermal expansion; always positive except below ^°C 

Kinematic eddy viscosity 

Thermometric eddy diffusivity 

Angular velocity of the vertical component of earth ^s rotation in 
Lake Michigan 

Width of the lake^ reference horizontal length dimension 

Depth of the lake, vertical dimension 

Total horizontal temperature difference between the edge and the 
central portion 
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TABLE 5-1 (Continued) 



Symbol Definition 



^ Stream function for the velocities 

P Fluid density at temperature T 

7 Depth to width ratio, 7 = D/L also nondimensional surface depth 

\ Tangent of the slope angle with the vertical X = tan (see Figure 

5-5) 

C Dimension of velocity C = (XgDAT/2QL 

p Nondimensional thermal Rossby number, p = agDAT/(2n)^ L^ 

G Nondimensional Ekman number, e = A/2nL^ 

Eddy Prandtl number, a = A/Aq 

MATHEMATICAL FORMULATION 

With all the considerations in the previous paragraphs, let us assume 
that the fluid is constrained in a long symmetrical, trapezoidal, cross- section 
model of the lake. It is also subjected to two oppositely imposed horizontal 
temperature gradients toward the center from both sides of the lake. The lake 
center is taken as the center of the Cartesian coordinates (Figure 5-2). The 
lake is rotating with angular velocity Q, (component of the earth's rotation nor- 
mal to the lake surface) about its vertical axis, which is anti-parallel to the 
gravitational force. The horizontal surface of the lake is considered free 
with no wind stresses. The mean flow is assumed to be zonal, thus the basic 
geostrophic thermal gradient relation is satisfied. All specifications of nomen- 
clatures are summarized in Table. 5-l» 

The governing equations of the lake model are then as follows: 

q^.V^q* + 2nk X q* + i V^p» - ag(T^ - T )k = AV»^J% (l) 

AqV^^T' - q*-V'T^ = 0, (2) 

V»-q^ = 0, (5) 



6h 



which are respectively the Navier-Stokes equations of motion, the energy equa- 
tion, and the continuity equation of the flow field. Note that the V and V*^ 
are two-dimensional gradient and Laplace operators. The above equations are, 
by axisymmetry, independent of Y. The boundary conditions are 



T' = T , q' = at z' = 0, (k) 

T' = -(1 - cos 2jtX')AT + T , ~^ = o at z = 7 (5) 

2 o oz 

||i- = , q' - at x= +1 L, (6) 

which imply that a horizontal temperature profile of the form of an inverse 
cosine curve is impressed upon the surface of the fluid of the lake, and the 
lake is insulated at both edges with only a slight conduction between the lower 
portion of the lake sides to the deeper layer of the earth. These differential 
equations together with all boundary conditions are coupled high order and 
nonlinear. Analytical direct mathematical solutions are very difficult. There- 
fore, some approximate approaches capable of making the equations more tractible 
are preferred. 



Nondimensional Equations 

It is convenient in the approximate method to introduce nondimensional 
variables and parameters. The following unprimed nondimensional variables are 
defined: 

rri t _ rp 

L ^ D ^ AT ^ ^ C C^p ^ 

o 

where C has the dimension of velocity. Since in our basic assumption the mean 
flow field is in geostrophic thermal gradient equilibrium, then the geostrophic 
thermal gradient relation: 

20 — = ag — 
dz "=" dx 

leads to 

agD 
c = 



^^-^'Sx _ ggPAT 



20, 2nL ' 

where D is the depth of the fluid. 
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Substitute all nondimensional variables into equations (l), (2), and (5) 
with the change of vertical dimension D to the same dimension as the horizontal^ 
L, but using 2=7, the depth to length ratio^ to denote the free surface, the 
governing equations become: 

p(q-V)q + k X q - kT + pVp = eS^q, (T) 



gV^T - crpq-VT = 0, (8) 

V-q = 0, (9) 

where 



p _ ^ ^ is the nondimensional thermal Rossby number as the (h) 



^ ^ of Hide^s (ISSk) parameter 



e = p is the nondimensional Ekman number 

a = A/Aq is the eddy Prandtl number 

Both p and e are very small number and a has the order of unity. The fluid is 
now bounded in 

--<x<7: ^0<z<7 where 7=7 , 



and boundary conditions are 

q - , T = at z = 0, (lO) 

■^-0 ,T=-(1- cos 2tix) at z = 7 (H) 

oz 2 

rVr 1 

q = on the side boundary and ■T: = Oatx = +— . (12) 

Since p is small, the ordinary perturbation method can be used. We may intro- 
duce a stream function by letting 
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t" = - JX Vt(x-Z) + €"M y (13) 



where the horizontal zonal velocity v is much greater than the meridional 
velocity, u, and the vertical velocity, w, which will appear in the following 
sections. 



Expand all the variables as powers of p as 



vn n (n) 

^ n=0 "^ ^ 



n=0 ^ 



n=0 ^ 



where the superscript n denotes the nth order of p- expansion. Equate terms in 
the powers of p. The zeroth p- order equations are: 

W°^ +v^°^ -T^°^ =0 (lU) 

Z X 

eV-;°) - A°) =0, (15) 

W°^ = 0. (16) 



The first p- order equations are 

^4(1) (1) (1) r (O) (O) (O) (O) (O) (O) (O) (O), 
Z X Z XZZ Z XXX X zzz X xxz 

(17) 

eW^) -^^'^ - [/°)v(°^ -^(°V(°)] = 0, (18) 

Z Z X Z Z 

eW") - CT[^^°^T^°) - ^^^°^t(°)] = 0, (IV) 

Z X X z 

and similarly for the other high p- order terms. 
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Boundary conditions are that all p- orders of velocities, and in turn, all 
p-orders of stream functions, should separately satisfy equations (lO), (ll), 
and (12); and that all p-orders of temperature should satisfy (lO) and (12) 
with (ll) for the zeroth p- order and 

T^^ = 0, (20) 



for other higher ith p-order temperatures, i = 1,2,3, etc. 
SOLUTIONS OF THE ZEROTH p- ORDER APPROXIMATION 

Temperature 

The zeroth p-order temperature can be obtained from equation (l6), which 
is simply the Laplace equation. The solution of equation (l6) satisfying boundary 
conditions (lO), (ll), and (12) is 

(o) Irz sinh 2nz ^ - /^ v 

r ^ = -L- - . ^ ^ cos 2TtxJ . (21) 

2 7 smh 2tC7 

The zeroth p-order temperature is purely conduction. 

Boundary Layer Scaling 

Equations (ih) and (15) are coupled and show boundary-layer characteristics 
because e is a very small number. In boundary- layer analysis the boundary- layer 
thickness should first be determined as a function of the small parameter, e. 
Equations (1^) and (15) imply that 

eW°^ + tl°^ = . (22) 



On the free surface and the rigid bottom we may scale the normal coordinates 
near the boundaries as 



cp :. €"^(z - 9) (25) 



where = on the bottom and = 7 on the surface. 
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Equation (25) substitutes into equation (22) with no stretching of the 
horizontal coordinate. The small parameter €, multiplying the highly differ- 
entiated term is important only near the boundary where the frictional force is 
the same order of magnitude as the other terms. This implies a = l/2. 

On the inclined side boundaries ( see Carrier 1955) 

where X is the slope of the slanted boundaries, and X = tan as shown in 
Figure 5-2; let 

I = ^""^(x + 0), (25) 

and 



where = +[2 - X(7 - z)] on each of the two side boundaries. By substituting 
equation (25) into equation (22), to the lowest order of (e), we obtain b = l/2. 
Therefore, all the boundaries have the boundary layer thickness of the order 

Following the procedure of boundary- layer additive type analysis, let 

(o) (o) (o) (o) (o) (o) ... 

r ^ = -^r ^ + ir ^ + ^r ^ + sr ^ + ^r . (26) 

where the left subscript denotes the region. The jt is dominating in the 
interior region and -t^^^, i = 1,2,5^^ ^^^ significant only on their respective 
boundaries. The various boundary- layer regions and their widths are shown in 
Figure 5-5' Coordinates for the regions of the lake cross- section are summarized 
in Table 5-2. 



Interior Solutions 

In the interior region, boundary effects are small. Equations (ik) and 
(15), with the first terms neglected, become: 

v(o) Jo) ^ ,(0) ^ 
^ ^ ^ = T^ ^ and ^\|;^ ^ = 
I X I z 
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JJL 

(I) 
12Y 




(I) indicates the Interior region 

(1) Indicates the surface region 

(2) Indicates the bottom region 

(3) Indicates the "eastern" boundary region 

(4) Indicates the "western" boundary region 
Figure 5-5. Boundary regimes of the Lake Michigan model, 



TABLE 5-2. Local coordinates in different regimes. 
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which yield: 



and 



v(o) l.cosh 27tz . >. n «/ \ 
E 2 smh 27t7 -^ ^ ^ ^ 



V°^ =B(x), 



(27) 



(28) 



where A(x), B(x) are integration constants and possibly functions of x. A(x) 
and B(x) can only be determined by matching up the interior solutions to the 
solutions obtained from the boundary- layers. 
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Boundary- Layer Solutions 

A. Region 1 — the upper boundary (free surface) 

In the upper boundary region, the stream function and the zonal velocity 
are, 

f - -j-f + i¥ ^ V - V + iv 

and the vertical coordinate is stretched as shown in Table 5-2. Equation (22), 
using the local coordinates, yields 

Since it^ ^ vanishes away from the upper boundary, no nonzero forcing term can 
be satisfied after integration of equation (29). Therefore equation (29) can be 
simplified as: 

Together with the appropriate boundary conditions, as 

5.^°^ ,(o) ^i^°^ , , 



or equivalently, 

equation (30) can be solved and we obtain: 

(o) cpi/f2 . , I- , , 

i\i/^ '^^ = -GTce^-^'^ sin 2Ttx cos cpi/v2 , (32) 



and 

( o) 1/2 Tte . ^ / cpi cpix 

iv^ = -eV^ ^^ sin 2^.x(cos ^ + sin ^) 
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(35) 



B. Region 2— the lower boundary (rigid bottom) 

Stretching the vertical coordinate^ equation (22) leads to an equation 
similar to (50). With the lower boundary conditions 

(o) (o) S^l;^''^ 

"^ ^""t = at cps = , (5U) 

we obtain 

2t' =-€:te ^^'^ sm 2Trx(cos ^ + sin ^) , (55) 



and 



M°^ = cV^vfl- «e-^^/^ sin 2«x cos ^ . 



(56) 



Matching the boundary solutions with that of the interior, the boundary 
conditions in equations (5I) and (^k) also yield the two integration constants 
in equations (27) and (28) as 

A(x) = sin 2Ttx(€V2 fz - - — ,^ ^ ) , 
^ ^ ' 2 sinh 2Tr7^ ^ 

and 

B(x) = -€tt sin 2Ttx . 



Then the interior solutions are: 

(o) 

j^ = Gjr sin 2:n:x , (57) 

(o) , 1/2 nr cosh 2jtz - 1 , . ^ ... 

_v' = (g / \r2 +•- — . , ^ ) sin 2nx . (58) 

I ^2 smh 21(7 

C. Regions 5 and h — eastern and western rigid side boundaries 

Using the local stretched coordinates normal to the side boundary, equation 
(22), to the lower order of €, leads to 
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Oil 



in region 5» 

With the boundary conditions, 

So) 



( o) St , , . 

v^ ^ = ^-— = at li = (UO) 

051 



we obtain. 



,(°) e(l + A.^)k kii/y/i" . ^ ,1 ,, ... cosh 2Ttnj - 1 r^n/ Mi • Hi\ 



(o) k|i/>/2' k|i . ^ a , . ,,cosh2itTii- 1 r?r 1/2-, 



where k = \ /(I + X ) . By symmetry, similar solutions are obtained in region 
k as 



4' 



(o) e(l + \^)k -)^iz/{2 . ^ a ,, , , , cosh2rtT]2-l r~,. kfo . k|p. 



,(o)_ -k|2/V2"_„ k 



4V • = -e --' ■ cos 



(111.) 



All the boundary- layer solutions are significant only on their respective 
local boundaries and go to zero away from their own boundaries. From equation 
(26), the general solution of the zeroth p- order approximation is: 
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.^^°^ = e[n sin 2^[-l . e^'^/^^ cos ^^ . e"^/'^ (cos -;#= . sin ^) ] 

_-.[..l/2-X(.-.)„ ^ (^^, 

J°) • o • r cosh 2itz - 1 ^ 1/2 ,r- rt -J^e^ '^^ rrr -z/\f2e z , 
^ = ^^^ "'^f 2 sinh 2.r """^ ^^W^^ -^^^ °°^-;|??] 

^ {sin 2n[i - A.(7 - ,)](<^°^^2^^- 1 + ^1/2^), ^^-k[x4-l/2-X(r-z)]/v[2r 

2 2 sinh 2tC7 

- cos kfx + 1/2 - X(y - z) 1 k[x-l/2+\(r-z) ]/\[2r kfx - l/2 + A.(y - z) 1 
cos ^ - e cos ^ ) ' 

and equation (21) is the zeroth p-order temperature. 



RESULTS OF THE ZEROTH p- ORDER SOLUTION 

For a small thermal Rossby number^ p^ and a high Taylor number^ equivalent 
to € ^ our solutions are within the range of lower symmetrical regime in a 
horizontal differentially heated rotating fluid (Fowlis and Hide I965). The 
governing equations are linearized by ordinary, perturbation expansion with p. 
Solutions are obtained by further singular perturbation expansion for each 
order of p in performing boundary- layer analysis with respect to the Ekman 
number €. 



Temperature 

The zeroth p-order temperature in equation (21) is obtained from heat con- 
duction, neglecting the convectional and the frictional effects. The temperature 
distribution is a linear function of depth peaked on the surface edge with some 
horizontal variation minimized at the central portion of the lake. Figure 5-^4- 
presents the horizontal temperature distributions of available data at the Great 
Lakes Research Division from 1965 to 1968 during the warm-up season. Figure 
5-5 shows the vertical temperature distributions during the same period of time 
in the central portion of Lake Michigan. Generally speaking, both the horizontal 
and the vertical temperature distributions show the same pattern as the solution 
indicated. The solution is in good agreement with the mean synoptic temperature 
measurement, especially in the interior region. On the boundaries, either hori- 
zontal or vertical, the agreement is not as good as the interior. On the upper 
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horizontal boundary, the real temperature decrease with depth is more rapid than 
the relation predicted by the solution. On the lower bottom layer, the real 
temperature decrease is much less. On the two slant side boundaries, real 
temperature gradients exist quite close to the shore, and in the central portion 
of the lake the temperature change is very small. The temperature variation is 
much greater than the solution yields. All these observations may indicate: 
(1) the impressed surface temperature is not an exact inverse cosine profile, 
and (2) the temperature is modified by convection as a result of the coupling 
of the temperature with the inertial terms in the equation of motion, particularly 
on the boundary layers. For mathematical simplicity, by smoothing the tempera- 
ture gradient, the inverse cosine temperature distribution is not too far from 
reality. As for the importance of convection in determining the temperature, 
higher p-order approximations should be attempted. Since the zeroth .p-order 
velocities have boundary- layer characteristics, and these velocities together 
with the zeroth p-order temperature gradients are forcing the first p-order 
temperature, it is expected that the higher p-order temperature will show 
boundary-layer characteristics and modification of the overall temperature dis- 
tributions. 



Velocity 

From the basic realistic assumption of the geostrophic thermal gradient 
balance, neglecting variable surface wind stresses, the zonal velocity dominates 
the overall circulation of the whole lake as expected. In the interior region, 
the pressure force balances the Coriolis force due to the zonal velocity together 
with the thermally induced body force. This basic assumption of geostrophic 
thermal gradient relationship (the balance of the vertical zonal velocity gradient 
and the horizontal temperature gradient) must be maintained. The interior solu- 
tion of the zonal velocity is shown in equation (38). It has the maximum magni- 
tude around one-fourth the width of the lake from the shore edge near the sur- 
face and approaches minimum as x approaches zero and half indicating the middle 
portion and both edges of the lake. During this warm-up season, on the eastern 
portion of the lake where x is positive, there is a northward zonal current. 
On the western portion of the lake, where x is negative in the horizontal coor- 
dinate of the lake model, there is a southward current. Both the horizontal and 
the vertical velocities are much smaller (in the order of e or ei/s) than the 
zonal velocity. From the interior stream function equation (57), there is no 
horizontal velocity (u = SVSz) in the interior region. The vertical velocity, 
w = - d\|f/Sx, leads to a broad but very small vertical circulation in magnitude 
in the interior region given as 

(o) P 

^w' ' = -€ 2Tr cos 2Trx ( J4.7) 

with a relatively strong sinking in the middle of the lake, and slight updwell- 
ing along the edges. 
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In region 1, that is^ in the free surface layer^ the zonal velocity is 
strengthened by adding to the interior zonal velocity a boundary zonal velocity 
of order (g-*-/^) in magnitude. The direction of the general zonal flow is the 
same as that of the interior, which shows a cyclonic circulation. The 
meridional velocity in the upper boundary is different from the interior where 
no east-west horizontal component develops. On the surface layer^ there is a 
meridional velocity 



..(o) St -1/2 ^ 1/2 ^ ^i/{^ . ^ / 9i . ^i\ 

^^' = -^- = € /^ -L_ = e/^-:^ e^^/^ sin 2^(cos ^ - sm ^) 



(^8) 



which flows from shore boundaries due east and due west^ toward the central 
portion of the lake. The vertical velocity in region 1 is much smaller than 
the horizontal velocity as the continuity equation implies. The vertical velo- 
city is: 

(o) (o) (o) (o) ^2 ^1/^2 cpi ^ .,^. 

W^ ^ = W^* + iW^ I^ ^ + G 2Tr e^^^"* cos ^ cos 27tK (1+9) 

which is the continuation of the interior velocity with only a slight increase 
in intensity. 

In region 2, the bottom layer, there is a meridional flow of an equal or- 
der of magnitude but opposite in direction to that of the surface, with a ver- 
tical velocity similar to that of the surface. The zonal velocity in the bottom 
layer is much smaller, as expected. It can be noticed that the residual zonal 
velocity near the bottom is the difference of the interior and the boundary 
zonal velocities. 

In the side boundaries, regions 5 sind U, the local coordinates (|^r|) are 
stretched normal to the inclined surfaces. From equation (^l) the meridional 
velocity in region 5, to the order of e, is: 

/ -1/2 . S S V ,( o) 





,(0) _ ^t^°^ 


that is 




(0) 


-eV^ k^(i + x2)^k|i. 



. k^i/>/2' . k^i . ^ , v.cosh 21^111- ^ rT\ 

)e sxn -f sxn 2^7 " ri3)( , ^^^^\^ - <^) 



(50) 



&k 



This meridional current flows away from the shore with a moderately high speed 
in the side boundary as |i is negative when X is positive, and vanishes in the 
interior region. By symmetry, a similar flow away from the west coast toward 
the center of the lake exists in region k. 

The local vertical velocity in the eastern alongshore region, to the order 
of e, is 



,w 



(o) 



ax 



= ._ei/2 k^d + x^) 



e^^^/ ^ sin^ sin 2TtX(r 



Tll)( 



cosh 2Trr]T - 1 
2 sinh 2ny 



(51) 



which represents sinking in the side bounda ry- layer . The local vertical velo- 
city in the alongshore boundary region is opposite to the vertical velocity in 
the interior. Therefore in the alongshore boundary region^ the overall vertical 
velocity may be small. Similar symmetrical phenomena exist in region k. 

The local zonal velocity in regions 5 and k, as shown in equations (^2) 
and {^h), are opposite in direction to the general pattern of the zonal circu- 
lation extended to the boundary from the interior. In the east and west along- 
shore boundary region^ the current direction is dependent on the difference in 
magnitude of the general zonal current and the local zonal current. Therefore^ 
the overall zonal currents in these two alongshore boundary regions may flow 
either southward or northward. Theoretically^ boundary currents very close to 
the shore should be smaller. 

In summing up the foregoing current analyses^ the general horizontal zonal 
circulation is shown in Figure 5-6. The meridional and the vertical circulations 
are shown in Figure 5-7a- When the overall temperature in the lake is less than 
i^.°C^. the coefficient of thermal expansion of water^ a, is negative. Hence the 
meridional and vertical circulation patterns of the lake are changed as approx- 
imately shown in Figure 5-Tt>* 
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Figure 5-6. General horizontal zonal currents. 
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Figure 5-7. Meridional and vertical circulation. 
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Since we have not analyzed in detail the magnitude of each component of 
velocity. Figures 5-6 and 5-7 only roughly show the possible general circula- 
tion pattern. A synoptic surface current computed from dynamic heights in 
June 1955, reproduced from Ayers et al. (1958), is shown in Figure 5-8. At 
the central portion of the lake, the general circulation pattern in Figure 5-6 
is in good agreement with that of Figure 5-8. A series of drogue studies for 
surface current measurement are shown in Figure 5-9 for the eastern shore area, 
and in Figure 5-10 for the western shore area. Tables 5-5 ^^^ 5-^ give the 
elapsed time for the drogue runs and the average drogue speeds for the east- 
shore and west-shore experiments, respectively. Generally, drogues drift toward 
the north near the eastern shore and toward the south near the western shore, 
which is in good agreement with the theoretical solutions. Winds are also in- 
dicated in Figures 5-9 and 5-10. In some cases, currents were flowing directly 
against the winds. In studying the thermal bar phenomena (the temperature gra- 
dient zone) in the Great Lakes, Rodgers (1965) indicated that the circulation 
associated with the thermal bar should be in the pattern as shown in Figure 5-11 
(Figure 7 of Rodgers* paper). When compared to our theoretical current pattern 
in a vertical section shown in Figure 5-7ti the major characters of these features 
agree well in general. Of course, our model does imply that this kind of circu- 
lation pattern exists only when the temperature in the thermal bar is lower than 
4°C. It is expected that when the sharp temperature gradient exists around the 
k°C zone, this kind of circulation pattern is more augmented by the acceleration 
of the vertical velocity in the boundary by the negatively buoyant force. 
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Figure 5-8. Surface currents in Lake Michigan, 29 June 1955 
(after Ayers et al. 1958)- 
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Figure 5-9 • Drogue studies for surface current mea,surements near 
eastern shore, (a) 2 May I968 (wind N, 10 kts?). (b) 2k June 
1968 (wind SW 15 kts). 
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Date 



TABLE 5-3 Drogue Studies on the East Side. 



Wind Drogue Station Start Finish cm/sec. 



2 May '68 N-10 



24 June 



SW-15 



A. 
B. 
C. 
D. 
E. 
F. 
G. 
H. 
I. 
J. 
K. 
L. 



0950 
1003 
1009 
1015 
1020 
1026 
1031 
1032 
1027 
1022 
1017 
1012 



1504 
1514 
1522 
1534 
1545 
1554 
1559 
1502 
1520 
1548 
1602 
1612 



7.14 
5.80 
5.80 
3.12 
6.24 
6.24 
4.01 
8.94 
5.80 
7.70 
10.26 
8.96 



30 July SE-20 



M. 


1007 


1623 


9.41 


N. 


1002 


1637 


9.41 


0. 


1005 


1328 


28.9 


P. 


0952 


1300 


19.5 


Q. 


Lost 






R. 


0931 


1220 


10.1 


S. 


0918 


1235 


12.7 



(c) 50 July 1968 (wind SxE, 20 kts?) 
Figure 5-9 (Concluded). 
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TABLE 5-4 Drogue Stations on the West Side. 



Wind Drogue station Start Finish cm/sec . 



10 May '68 Calm 



Miles 



27 June 
87^50' 
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Figure 5-10. Drogue studies for siirface current measurements near 
western shore, (a) 10 May I968 (wind slight), (b) 27 June I968 
(wind NxW, 6 kts). 
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Figure 5-11. Circulation associated with "Thermal Bar" 
(after Rodgers I965). 



THE FIRST p- ORDER SOLUTIONS 

As we have already mentioned, the first p- order governing equations of the 
model are equations (1?)^ (18)^ and (19). The forcing terms which are the 
products of the zeroth p-order temperature gradients and the velocities (as 
shown in equation (19))^ are of the order unity in the interior region, of 
the order (e"""/^) for the individual local boundary- layer regions at their re- 
spective boundaries, and with another forcing term of the order unity due to 
the overall interior influence. Since the forcing terms are of the boundary- 
layer type, it is expected that the first p- order temperature will have strong 
boundary- layer characteristics. In each of the boundary regions, the first p- 
order temperature can be separated into two parts; one part is defined by the 
overall existing zeroth order interior stream function (^t^^^), and the other 
part is due to the local stream function which is important only in its own 
boundary region. In this way the boundary- layer portion which is forced by 
the terms of order of (e ^/^) can be treated with boundary- layer analysis. 
After applying the respective local boundary conditions, solutions are obtained 
for all the first p- order temperatures influenced by local zeroth p- order stream 
functions. Then at all the new boundary conditions, solutions are obtained by 
applying the first p-order solution which has the boundary- layer characteristic 
as a perturbation of that part of the temperature under the Laplace operator 
(which has no boundary- layer characteristic). In the interior region the tem- 
perature is separated as a general solution and a particular solution, with the 
forcing terms of the order of unity operating only on the particular solution. 
After transferring the nonhomogeneous terms as a particular solution into the 
new boundary condition and matching with the boundary solutions to the order of 
unity, the overall first p-order temperature can be solved. The result of the 
first p-order temperature turns out to be of the order of unity and there are 
some modified variations of the order ( e^/^) in the boundary regions. 
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After obtaining the first p-order temperature, the first p-order stream 
function and the zonal velocity can be obtained from the coupled equations (I7) 
and (18) following the method used in the zeroth p-order solutions, and with 
more complexity. 



CONCLUSIONS 

After the winter mixing, water in Lake Michigan is almost homogeneous. 
The spring warming season begins during late April or early May, ' Beginning at 
this time the lake starts to warm from the shallow water along the shore 
progressively outward to the main body of water which is still at k^'C or lower. 
During this period, the land and air are much warmer than the water. A stable 
situation usually exists in the atmosphere above the lake. It is expected that 
much less energy is transferred from the wind to the lake water. A sharp tem- 
perature gradient always exists within a few miles along the shore which leads 
to a balanced geostrophic thermal gradient model of Lake Michigan in a simpli- 
fied long trapezoidal shape. A double perturbation method is used in obtaining 
approximate solutions from these coupled nonlinear momentum and energy equations 
by expanding first with the small thermal Rossby number, p, and then, to each 
order of p- using singular perturbation, with respect to the small Ekman number, 
€. All boundaries are found to have the same order of boundary- layer thickness 
in the order of ( € ^f). The zeroth p-order temperature is found to be due en- 
tirely to conduction. The first p-order temperature shows that convection is 
important and it has boundary- layer variations. All the velocities do have 
boundary- layer characteristics. The general zeroth p-order horizontal zonal 
velocity shows a cyclonic circulation pattern with variations in the along- 
shore boundaries. When the overall temperature of the lake is still lower than 
k°C, the horizontal meridional velocities and the vertical velocities show a 
broad upwelling of a small order of magnitude in the middle portion of the lake 
and a sinking zone along the edge of alongshore boundary- layers coupled with the 
divergent horizontal surface meridional current and the convergent bottom cur- 
rent. This forms two gigantic but weak counter-rotating cellular motions. Two 
horizontal meridional currents flowing outward from their respective shores in 
the two alongshore boundary regions to joint the thermal bar convergences sug- 
gest the possible existence of another upwelling very close to the shore bot- 
toms. Though detailed characteristics have not been wholly demonstrated, the 
general circulation pattern compared with results obtained by Ayers et al. 
(1958) sind Rodgers (1965) indicate very good agreement. From a physical point 
of view, in the spring warming season the shallow water alongshore region re- 
ceives a great amount of heat energy from the air and the land. The temperature 
in these regions increases very rapidly while that of the main body of water 
does not respond as quickly and remains at the temperature of maximum density. 
This temperature gradient induces an inclined surface water slope tilted up 
toward the land, which induces a "Hadley cell" (1755)- The thermal current is 
originally flowing outward from shore but is deflected. by the Coriolis force 
which results the northward or southward zonal current. Due to the relatively 
small dimensions, the boundary effects are more important in lakes than in 
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oceans. Though we have not analyzed the numerical magnitude of all variables^ 
the general circulation pattern has demonstrated all the main expected results, 



FUTURE PIAN 

Our preliminary model is too rough to obtain detailed characteristics of 
the lake circulation either in magnitude or in flow pattern. It is my hope to 
be able to further refine the model both in dimension and magnitude of some var- 
iables in similation of the prototype of the lake. For example^ the horizontal 
dimension and the vertical dimension should be scaled separately; the vertical 
velocity should have a different velocity scale than the horizontal one; the 
horizontal eddy viscosity and the vertical eddy viscosity should have different 
mean values. As in perturbation procedures, higher orders of approximation are 
very important. The detailed analyses of the first g-order solutions and the 
estimation of approximate error from the second P-order equations are scheduled 
to be carried out. The combined final solutions should be programmed into com- 
puter calculations and the results should be compared with all available field 
data including some of the current meter buoy stations of the Great Lakes- 
Illinois River Basins project in 196^. Furthermore, I also intend to consider 
the thermal body force as a slowly time dependent function and to extend the 
application of the model for the entire year. 
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